Betti number signatures of homogeneous Poisson point processes. 
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The Betti numbers are fundamental topological quantities that describe the fc-dimensional connec- 
tivity of an object: [3o is the number of connected components and f3u effectively counts the number 
of k- dimensional holes. Although they are appealing natural descriptors of shape, the higher-order 
Betti numbers are more difficult to compute than other measures and so have not previously been 
studied per se in the context of stochastic geometry or statistical physics. 

As a mathematically tractable model, we consider the expected Betti numbers per unit volume 
of Poisson-centred spheres with radius a. We present results from simulations and derive analytic 
expressions for the low intensity, small radius limits of Betti numbers in one, two, and three di- 
mensions. The algorithms and analysis depend on alpha-shapes, a construction from computational 
geometry that deserves to be more widely known in the physics community. 

PACS numbers: 02.50.Ey Stochastic processes; 02. 40. Re Algebraic topology; 05.10.Ln Monte Carlo studies 
Keywords: Topological invariants, Betti numbers, Euler characteristic, Poisson-Boolean model, alpha-shapes, 
continuum percolation. 



I. INTRODUCTION 

Topological measures of shape are finding increasing 
use in the study of point or coverage processes and the 
characterisation of complex three-dimensional structures 
0. This is because topology is independent of geome- 
try, and so both types of information are necessary to 
fully characterise spatial structure 0, |j| . The most com- 
monly studied topological invariants are the number of 
connected components (the zeroth order Betti number, 
Po) and the Euler characteristic (\, the zero-dimensional 
Minkowski measure from integral geometry). This pa- 
per also investigates Pi and the higher-order Betti 
numbers that count the number of independent handles 
(non-contractible loops) and enclosed voids. 

The Betti numbers are closely related to the Eu- 
ler characteristic via the Euler-Poincare formula: \ = 
Po — Pi + • • • — Pm- For subsets of K, the Euler char- 
acteristic is exactly the number of components. In K 2 , 
there are only two independent quantities from the three, 
since x — Po — Pi- Thus, as x an d Po are already well- 
known quantities in statistical physics, the higher-order 
Betti numbers give intrinsically new information only in 
dimensions three and higher. Nonetheless, it is instruc- 
tive to study the Betti numbers directly in both two- and 
three-dimensions as they give a more direct description of 
the topology than the Euler characteristic. For example, 
recent work on the 2D Griffiths' model has used Betti 
numbers of the different states to characterise the phase 
transition p|. 

In this paper, we present analysis and simulations of 
the Betti number signatures of Poisson point patterns. 
The Poisson point process is the most widely studied 
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model in stochastic geometry and is frequently used as a 
null hypothesis for comparison with physical systems [j| . 
In general, a signature for a point pattern is defined by 
attaching spheres of radius a to each point and comput- 
ing some quantity of interest as a function of a. Thus, 
the Betti number signatures contain both topological and 
geometric information about the distribution of points in 
space. In applications, such signature functions can be 
used to detect differences between simulations and physi- 
cal data, or to provide insight into the physical processes 
that generated a particular distribution of points. Ex- 
ample applications will be the topic of a future paper. 

We give a brief overview of the simulation of Pois- 
son point processes in Section III Al Our computation 
and analysis of the Betti number signatures use alpha- 
shapes [a, — a construction from computational ge- 
ometry that is dual to the union of spheres of radius a. 
The alpha-shape is a subcomplex of the Delaunay trian- 
gulation of a set of points, so we can draw on extensive 
results about Delaunay complexes of Poisson-distributcd 
points. We summarise the alpha-shape and Betti number 
algorithms in Sections III Bl and III CI Results of the sim- 
ulations are presented in Section II I II The final section 
(|IV|) of the paper gives derivations of the low-intensity 
small-radius behaviour of the Betti numbers of Poisson- 
distributed spheres. 



II. SIMULATION METHODS 

A. Poisson point processes 

A Poisson point process in M. d with constant intensity 
A is easily simulated in the unit c£-cube by generating N 
points with d coordinates chosen from a uniform random 
distribution on [0,1]. The number of points, N, is a 
random variable generated from a Poisson distribution 
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with mean A, 

Pr(N = n)= X n e- x /nl 

For large values of A, Poisson distributed numbers are 
well approximated by a normal (Gaussian) distribution 
with mean A and standard deviation VA, i-e., 

/•n+0.5 

Pr(N = n) = / f{x)dx, 

Jn-0.5 

where f(x) is the normal probabiltiy density function 
/(*) = -2=exp[-(x-A) 2 /2A]. 

V ZlTA 

Given a realisation of a Poisson point process in the 
unit e?-cube, label the N points Xi, X2, . ■ ■ , Xn and 
place identical balls of radius a centred at each point, 
Bd{Xi,a). The Betti number signatures are defined to 
be 

Pk{a) = fa (\J B d (X u a)^J for k = 0, 1, . . . , d - 1. 

Algorithms for computing the Betti numbers are de- 
scribed in the following section. Expected values per unit 
volume, Ef3k(a), are estimated as mean values calculated 
from many independent realisations of points in the unit 
d-cube. 



B. Betti numbers of alpha shapes 

The union of balls of radius a centred at the points 
X\ , . . . , Xn has a geometric dual called the alpha- 
shape that is a subset of the Delaunay triangulation of 
Xi, . . . , Xjy. The nerve theorem of topology guarantees 
that the Betti numbers of the union of balls are identical 
to those of the dual alpha-shape 8J. Since the alpha- 
shape is a discrete simplicial complex, the Betti numbers 
are computable via linear algebra techniques for data in 
any dimensions 9| • The classical algorithm is impractical 
for large complexes however, and for points in one, two, 
or three dimensions there are more effective geometric 
approaches. 

In ID there is only the number of connected compo- 
nents to consider, Po(a), and this is determined entirely 
by distances between adjacent points. For points in 2D 
and 3D, the Betti numbers may be computed via an 
incremental algorithm due to Delfinado and Edelsbrun- 
ner ^(j that gives Pk(ct) at all values of a. The essential 
aspects of their approach are as follows. 

Firstly, the simplices of the Delaunay complex — the 
vertices, edges, triangles, and so on — are ordered by 
the radius of the smallest sphere that touches the points 
of the given simplex and contains no other data points. 
This radius is called the alpha-threshold, ax- If more 
than one simplex has the same alpha-threshold, they are 



ordered from lowest dimension to highest. The ordering 
of simplices {a\, 0%, . . . , <j n } such that otT{o~i) < a r(c :) ) 
if i < j is a filtration. A sequence of subcomplexes Cj — 
Ui=i °i i s now built by adding one simplex at a time. 
Each fc-dimensional simplex either creates a new fc-cycle, 
or destroys a (fc-l)-cycle. For example, when an edge is 
added, it either generates a loop or connects two disjoint 
components. The Betti numbers of Cj+i are related to 
those of Cj by: 

Pk(j + 1) = Pk(j) + 1 if creates a /c-cycle 
Pk-i(j + 1) = Pk-i(j) - 1 if ctj+1 destroys a (fc-l)-cycle. 

The problem of determining whether a fc-simplex cre- 
ates a fc-cycle is non-trivial in arbitrary dimension d. 
Fast algorithms based on union-find data structures are 
possible when k = 1 and, by Alexander duality 0, 
k = d — 1. Thus, this incremental approach is effective 
only for points in d — 2, 3 dimensions. The duality argu- 
ment requires the complex to be a subset of the G?-sphere, 
but this is easily accounted for by adding a point at in- 
finity to the Delaunay complex, and an extra fc-simplex 
for each (fc-l)-face on the convex hull. 

Each fc-simplex in the filtration, as it is added to the 
complex, is marked +1 if it is found to create a fc-cycle, 
and —1 if it destroys a (fc-l)-cycle. Then the Betti num- 
bers of the alpha-shapes are calculated as 

Pa(a) = #{+1 vertices < a} — #{ — 1 edges < a}, (1) 
Pi(a) = #{+1 edges < a} — #{ — 1 triangles < a}, 
/32(a) = #{+1 triangles < a} — #{—1 tetrahedra < a}. 

Note that these signature functions may be evaluated at 
arbitrary fineness in a, with no additional computational 
cost or complexity. A single traversal of the marked sim- 
plices, in the filtration order is all that is required. 

C. Periodic boundary conditions 

To avoid boundary effects from restricting the domain 
to the unit d-cube, we build the Delaunay complex with 
periodic boundary conditions. 

For points in 2D, this means the Delaunay complex is 
a triangulation of the 2-torus, and the Betti numbers for 
sufficiently large a are 0q — l,/3i = 2, /3a = 1. The in- 
cremental algorithm for computing the Betti numbers is 
still valid, provided the final triangle added to the filtra- 
tion is marked as creating a 2-cycle (which we know it 
must, a priori). When a is below the percolation thresh- 
old, /?i(a) may be interpreted as the number of holes per 
unit area. If we use periodic boundary conditions, then 
above the percolation threshold 0i(a) includes the cy- 
cles around each axis of the torus, and so the number of 
holes in the unit square, as one would intuitively define 
them, is really (3\(a) — 2. This issue is related to the 
problem of whether or not to count the spanning cluster 
when studying the connected components in percolation 
theory. 
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For points in 3D, periodic boundary conditions invali- 
date the algorithm for the determination of 2-cycles. The 
topology of a unit cube with opposite faces identified is 
that of a 3-torus: /3q = = 3, 02 = 3, /?3 = 1. Thus, 
the sequence of subcomplexes obtained from the filtra- 
tion are subspaces of the 3-torus, not the 3-sphere, and 
the Alexander duality theorem no longer applies. 

In practice, we apply the duality algorithm to the 2- 
cycle detection problem anyway. The result is that ex- 
actly three faces (triangles) from the filtration are in- 
correctly identified as "destroying 1-cycles" when in fact 
they create the three 2-cycles that are homologous to 
the three coordinate planes of the periodic cube. This 
means there are three extra triangles marked —1, and 
three fewer marked +1, than there would be if we had a 
direct algorithm for detecting 2-cycles. All other edges, 
faces, and tetrahedra are correctly marked ±1, provided 
the final tetrahedron is identified as creating a 3-cycle 
(which we again know a priori). Comparison with the 
formulas given in JQ) shows that for both k = 1,2, 

/3 c k alc (a) = Pl rue (a) - #{mislabelled triangles < a}. 

If we consider the alpha complex as a increases, it should 
be clear that the mislabelled triangles have the small- 
est alpha-thresholds for which each of the three torio- 
dal 2-cycles exist. This represents the second percolation 
threshold: a critical radius, above which, with prob- 
ability one, the unoccupied space no longer percolates. 
Thus, for the mean values of the first and second Betti 
numbers we have (for k = 1, 2) 

Ef3 c k alc {a) = Ef3l rue (a) for a < a 2 

Epi alc {a)=E(3l rue {a)-? l for a > a 2 . 



D. Implementation 

There are three publicly available implementations of 
alpha-shapes: Edelsbrunner's gro up fill . Clarkson's hull 
code ^3 > an d the CGAL library |l3t Il4[ . None of these 
has provision for periodic boundary conditions, and only 
the first has support for computing the Betti number 
signatures. The CGAL library (written in C++) has the 
most general interface, so we use the CGAL implemen- 
tation of two- and three-dimensional Delaunay triangu- 
lations and alpha shapes and extend it as follows. 

First, N points in the unit square or cube are gener- 
ated with uniform random coordinates, and each point 
is given a label. Periodic boundary conditions are simu- 
lated using translated copies of the original data points 
with each translated copy of a point given the same label 
as the original. The simplest approach to generating the 
translated points is to map all the original data points to 
the 8 adjacent squares in 2D, or the 26 adjacent cubes 
in 3D. This creates a significant overhead in the number 
of points to be triangulated — 9A^ and 27N respectively. 
There is also an increasing degree of redundancy for large 



N, since almost all of the translated points have no ef- 
fect on the triangulation within the original domain. For 
reasonable numbers of points (N > 50) we can there- 
fore reduce the overhead by translating only the data 
points in the appropriate half-cube along each axis, lead- 
ing to AN and 8A^ points to be triangulated in 2D and 
3D respectively. The minimum requirement on trans- 
lating points that guarantees a correct triangulation is 
in principal even less: only points that belong to a De- 
launay cell whose circumsphere intersects the boundary 
of the unit square or cube need to be translated to the 
opposite side [l^. However, we find that the increased 
complexity of this approach outweighs any saving from 
the reduced number of translated points. 

The second step is to build the Delaunay complex and 
alpha shape on the enlarged set of data points using the 
CGAL library routines. We must then identify the ele- 
ments of the Delaunay complex that comprise the peri- 
odic domain. The criterion we use is that the centroid 
of the cell (or face, or edge) is either interior to the unit 
cube, or lies on one of the x = 0, y = 0, or z = planes. 
The topological integrity of the Delaunay complex with 
the periodic boundary conditions is checked via the labels 
attached to the vertices. 

Finally, we implement a filtration data structure and 
the incremental Betti number algorithm as described in 
Section lll Bl The CGAL alpha shape data structure gives 
us direct access to the alpha-thresholds of each simplex 
(i.e., the cells, faces, and edges), so this is relatively 
straightforward. The C++ code is available from the 
author on request. 

The simulations reported in Section lTTTI were performed 
on a PC with Intel Pentium 4 processor. The two- 
dimensional simulations involved 1000 realisations with 
A = 10 5 and ran overnight. The three-dimensional sim- 
ulations involved 50 realisations with A = 10 5 and took 
five days. The dramatic increase in time for the 3D simu- 
lations is due to the intrinsic additional complexity of 3D 
Delaunay complexes and alpha shapes, the extra points 
needed to simulate periodic boundary conditions, and the 
need for both a forward and backward traversal of the fil- 
tration to mark the simplices. 



III. RESULTS AND ANALYSIS 

In this section, we summarise theoretical results and 
compare these to data obtained from computer simula- 
tions of Poisson point processes in two and three dimen- 
sions. 



A. ID 

The expected number of components per unit length 
in a ID Boolean model is well known ^(|. If the intensity 
of the Poisson-point process is A, and the shapes are line 
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FIG. 1: Expectations per unit d-volume of the Euler charac- 
teristic, x/A, as a function of the reduced density, r\ = ajd,a d \, 
for d — 1,2,3. The quantity u>d is the d-volume of the unit 
d-sphere. 
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FIG. 2: Deviation from theory of the computed mean Euler 
characteristic per unit d-volume as a function of the reduced 
density, 77, for d = 2, 3. In the 2D case, the mean Euler 
characteristic is computed from 1000 simulations of a Poisson 
point process with A = 10 5 in the unit square with periodic 
boundary conditions. In the 3D case, we have used 50 simu- 
lations with A = 10° in the unit cube with periodic boundary 
conditions. 



segments of length 2a then 

Ex(a) = E/3 (a) = Ae" 2aA . 

This result is included for completeness and ease of com- 
parison with the two- and three-dimensional cases, see 
Figure ^ 

B. 2D 

For discs of radius a centred at points from a 2D Pois- 
son point process of intensity A, we study the expecta- 
tion per unit area of the following topological quanti- 
ties: the number of components, /3o(a), the number of 
independent cycles, /3i(a), and the Euler characteristic, 
X = Po — Pi- For our simulations, we use an intensity 
of A = 10 5 in the unit square, and compute mean values 
of the Bctti numbers from 1000 realisations. Results are 
presented in Figures and 0] In the plots of these figures 
we mark the 2D continuum percolation threshold from 
[l7l | of rj c — 1.1280586. The critical value is included as 
a reference point only, since the Betti numbers are not 
sensitive indicators of percolation. 

The expectation per unit area of the Euler character- 
istic is known from stochastic geometry to be 0, 0] 



E x (a) = A(l - na 2 X)e-' KC 
= AOL-^e-". 



(2) 



This expression is more naturally a function of the re- 
duced density, 77 = 7ra 2 A, and we often use 77 as the in- 



dependent variable rather than the radius a. The differ- 
ences between the expression @ and the computed mean 
values of the Euler characteristic are less than 10~ 4 , and 
decrease as r\ increases, see Fig. |21 

The connected components of randomly distributed 
overlapping discs are studied extensively in percolation 
theory. It is common in this context to express the ex- 
pected total number of components per unit area as the 
sum 



k=0 

where pk is the expected number of k-mers per unit area 
(a fc-mer is a cluster built from k discs). Although there 
are no known analytic expressions for E(3q as a function of 
disc radius cv, integral expressions for pf. and low-density 
expansions are given in |18| . The expansions for pk are 
given in terms of the reduced density for the limit 77—5-0 
and presented in Tabled for reference. Their sum gives 



E(3 (r))/\ = 1 - 277+ 1.564177 2 



0.687877^ 

(3) 



+ 0.219777 4 + 0(rn. 

A comparison between this expansion and the computed 
mean values obtained from simulations is shown in Fig. El 
There is extremely close agreement for 77 < 0.5. 

An expansion for Efii may be deduced from the ex- 
pressions for Efto and Ex above. However, we make an 
independent analysis of the shape of Poisson-Delaunay 
cells in Section IIV Bl and find that for small 77 



EfoW/X = 0.0640r/ 2 + 0(t7 3 ). 



(4) 
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TABLE I: Coefficients in the expansions of pk for the 2D 
Poisson-Booiean modei of discs with radius a for the limit 
t] — Tv\a 2 — > 0. Results are from fl^l . 



n° v 1 



v 



if 



pi/X 1 -4 8 -10.6667 10.6667 
p 2 /A 2 -11.3079 32.2915 -62.0415 
p 3 /X 4.8720 -35.3346 129.6895 
p 4 /X 13.022 -114.823 
gs/A 36.728 
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FIG. 3: 2D Betti numbers. Results from 1000 simulations 
of, on average, 10 5 points in the unit square. Mean values 
of the Betti numbers per unit area, /3o (blue dots) and /3i 
(magenta dots) are given as functions of the reduced density, 
■q (main) and the disc area fraction = 1 — e~ v (inset). The 
percolation threshold is marked by the dotted vertical line at 
T] c = -1.1280586, or equivalently 4> c = 0.676339. 



Our simulation data show that this leading order be- 
haviour holds for i] < 0.3, see Fig. 0] 

The logarithmic axes used in Fig. 0] show that E(3q / A 
levels out at 10~ 5 . This is exactly as expected since for 
large radius, the alpha-shape has one connected compo- 
nent, and the value of A is 10 5 for these simulations. As 
discussed in Section III CI periodic boundary conditions 
mean that for sufficiently large radius we know 0i = 2. 
Thus we would expect to see Ef3i / A level out at 2* 10~ 5 in 
Fig. but the range in this plot does not extend to large 
enough r/. This shows that periodic boundary effects are 
negligible for the data from these simulations. 



3D 



In the three-dimensional Poisson-Booiean model of 
balls with radius a, the relevant topological quantities 
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FIG. 4: 2D Betti numbers. Exactly the same data as in 
Fig. [J] but plotted here with logarithmic axes to emphasise 
the quadratic scaling of /3i at small r\. The solid black line 
shows the leading order behaviour E/3i/X ~ 0.0640?; 2 derived 
in Section IPTBI 




FIG. 5: 2D connected components. Here we compare the 
computed mean values (blue dots) of the number of connected 
components, /3o/A with the theoretical expansion (pale blue 
line) for small rj given in @. 



are the number of components, Po(a), the number of in- 
dependent handles, (3\{a) 1 the number of enclosed voids, 
(3i(a), and the Euler characteristic, x = Pa — Pi + (3% ■ For 
the simulations we use an intensity of A = 10 5 in the unit 
cube and compute mean Betti numbers from 50 realisa- 
tions. Results are presented in Figures|Bland[7| We again 
mark the continuum percolation thresholds in these plots 
as reference points. Recall that in three-dimensional per- 
colation there are two critical densities: 771 = 0.341889 
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TABLE II: Coefficients in the expansions of pu for the 3D 
Poisson-Boolean model of balls with radius a for the limit 
77 = |-7rAa 3 — > 0. Results are from 18]. 



77° 77 1 



if 



pi/X 1 -8 32 -85.3333 170.6667 

p 2 /X 4 -49 302.2238 -1250.5030 

pz/X 22 -359.4203 2959.1209 

Pi/X 139.7867 -2842.60 

p 5 /X 964.68 



|i9j . is the point above which a span ning cluster exists 
with probability one, and 772 = 3.5032 [2(j is the density 
above which the unfilled space no longer percolates. 

The expectation per unit volume of the Euler charac- 
teristic is again known from stochastic geometry |5j to 
be: 



E X (v) = A(l - 377 + 



32 



(5) 



where r\ is the reduced density 77 = |7rAa 3 . Our com- 
puted mean values match this expression closely, with 
differences less than 10 ~ 3 and decreasing with 77 as shown 
in Fig. H 

As for the 2D model, the expected total number of 
components per unit volume may be expressed as the 
sum of numbers of fc-mers. Integral expressions and low- 
density expansions for the expected number of fc-mers per 
unit volume, pu, are given in an d repeated here in 
Table ITT1 From these expansions we find that for 77 — > 



E0o(r))/\ = 1 - 4ry + 5r) 2 - 2.7431 rf 

+ 1.3646 7] 4 + 0(ri 5 ). 



(6) 



The computed mean values match this expansion ex- 
tremely closely for r\ < 0.3, see Fig. |SJ 

The leading order behaviour for (3\ and 02 is derived 
from the Poisson-Delaunay analysis in Sections II V CI and 
II V Dl where we show that for small 77 



Epi(v)/X 



0.5747 rf 
0.015 t? 3 - 



f 0( ?7 3 ) 



(7) 
(8) 



Again, the computed mean values show exactly this lead- 
ing order behaviour for 77 < 0.3, see Fig.0 

Recall from Section III CI that with periodic boundary 
conditions and 77 > 772, the second percolation threshold, 
there is a systematic error in the computed mean values 
of (3i and j3^- For the data presented here, the error is 
3* 10~ 5 , which is three orders of magnitude less than the 
value of E(3\{r)2)/\ and four orders less than Efaiift)/^- 
A close inspection of Fig. [7| shows that this error is sig- 
nificant only for the computed values of Ef3\(rf)/\ with 
?7 > 6. 
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FIG. 6: 3D Betti numbers. Results from 50 simulations of, on 
average, 10 s points in the unit cube. Mean values of the Betti 
number per unit volume, /3o (blue dots), /3i (magenta dots), 
and P2 (red dots), are plotted as functions of the reduced 
density 77 = §7ra 3 (main) and ball volume fraction (f> = l — e~ v 
(inset). The two critical densities from percolation theory 
are marked by dotted black lines at 771 = 0.341889 (cj>i = 
0.289573) and r/ 2 = 3.5032 {<f> 2 = 0.9699). 
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FIG. 7: 3D Betti numbers. The same results as in Fig. El but 
plotted with logarithmic axes to show the power-law scaling of 
/3i and fii for small 77. The solid black lines show the leading 
order behaviour of EfhJ X ~ 0.5747 77 2 , and Ef3 2 /X ~ 0.015?? 3 
derived in Sections IIV CI and II V Dl 



IV. POISSON-DELAUNAY CELL ANALYSIS OF 
ALPHA SHAPES 

The probability distribution for the size and shape of 
a cell in the Delaunay complex of a Poisson point process 
is completely characterised by a result due to Miles [2l| 



7 




n 



FIG. 8: 3D connected components. The computed mean val- 
ues of Po/X (blue dots) compared to the low- intensity expan- 
sion (pale blue line) given in ©. 



on its boundary or in its interior. This is the circum- 
sphere of Xq , . . . , X m ; let c and r denote the circumcen- 
ter and circumradius respectively. The vertices of a PDC 
are given in vector form by X$ = c + riij , where Uj is a 
unit vector pointing from the centre to the point X{. The 
circumradius is a measure of the size of a simplex, and 
the unit vectors specify its shape. 

The distribution of PDCs is completely specified by the 
following probability density function (pdf). This result 
is due to Miles [21| , and implies that the circumradius of 
a PDC is independent of the positions of its vertices; 

h m (r, Uo, • . . , u m ) = o(A, m)A m r m2 ~ 1 exp(-Aw m r m ). 

(9) 

The constant cj m = it™ 1 / 2 /T{m/2 + 1) is the volume of 
the m-dimensional unit sphere, and 

_ 7r( m2+1 )/ 2 r(m 2 /2){2Ar[(m + l)/2]} m 
a( ,m '~ m™- 2 r(m/2) 2m + 1 r[(m 2 + l)/2] ' 



and given in Eq.©. The criteria for a simplex from the 
Poisson-Delaunay complex to belong to an alpha shape 
are based only on the size and shape of that simplex, 
and that of its adjacent simplices. The ergodicity of 
the Poisson-Delaunay complex means that the expected 
number of fc-dimensional simplices, a, in a bounded re- 
gion, R, that satisfy condition A, is related to the prob- 
ability that a randomly selected simplex has property A: 

E#{a G R | a is A} = X k \\R\\Pr{A), 

where X k is the intensity of the fc-dimensional cells, not 
the vertices (which have intensity Ao = A). Since the 
Betti numbers of alpha-shapes are determined by num- 
bers of simplices with certain properties, see ||TJ, the 
Poisson-Delaunay cell (PDC) distribution can be used 
to obtain results about the Betti numbers of an alpha- 
shape. 

This section summarises the relevant results about the 
PDC distributions in two- and three-dimensions, and 
then derives low-intensity expansions for the expectation 
per unit area of f3\ in 2D and expectation per unit volume 
of Pi and fa m 3D. 



A. Distributional properties of PDCs 

For an extensive review of Poisson Delaunay cells, see 
the book Spatial Tesselations (2nd ed., Section 5.11)p2|. 

We start by considering a Poisson point process with 
intensity A in R m . A Poisson Delaunay cell is an Tri- 
dimensional simplex, i.e., the convex hull of to + 1 
points Ao, . . . , A m from the Poisson point process, such 
that there exists an (m-l)-sphere that has each point 
Aq, . . . , X m on its boundary and no other points either 



The dependence of h m on the is hidden in the func- 
tion A m , defined as the volume of the m-simplex with 
vertices at uo, . . . , u m . It is therefore a constant, 1/(to!), 
times the determinant of a square matrix with m rows 
containing the vectors (iij — Uo), for i = 1, . . . , to. 

Various distributional properties of PDCs can be de- 
rived from this pdf. In particular, Muchc [HHil has sim- 
plified the pdf for the two- and three-dimensional cases, 
finding in 2D: 

/a(r, 0i , 2 )=2(^A)V exp(-A7rr 2 ) • 

2 . 01 . 02 . 01+02 ( 10 ) 



where < r < oo is the circumradius, < 2 < 2ir — 
0i, and < 0i < 2ir are the central angles AocAi and 
AicA 2 . 

In 3D, we can choose a coordinate system so that 
the circumcentre is at the origin, and three points 
(Ai, A 2 , A3) of the tetrahedron lie in the plane x = cos 6, 
where 8 is the angle between the normal to this face (i.e., 
the positive x-axis) and one of its vertices. The y- and z- 
coordinates of the vertices in this triangular face are then 
determined by the central angles 0i and 02- The distri- 
butional properties of this face are those of a "typical" 
face in a Poisson-Delaunay complex. The fourth vertex 
of the tetrahedron is specified by the height of the tetra- 
hedron, h, and an angle 7. Muche |23| showed that the 
pdf for a Delaunay tetrahedron separates into factors: 

/ 3 (r, 6, h, 0i, 02, 7) = f R (r)fe,H(9, h)M<fn, h)Ml) 

(11) 
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with marginal densities: 
32^ 3 A 3 



/n(r) 

Mi) 



o 47rA o. 
r exp( — r J, < r < oo, 

64 

0</i<l + cos6>, < 9 < tt, 



16 / . 01 . 02 . 01 + 02 

— =r sin — sin — sin 

3tt 2 V 2 2 2 



O<0 2 <2tt-0i, < 01 < 2tt, 
1 

27' 



< 7 < 2tt. 



B. Empty triangles in 2D 



We now derive conditions on a Poisson Delaunay cell 
in M 2 that guarantee the 2-simplex is excluded from the 
alpha-shape, but all its edges are included in the alpha- 
shape. This implies the existence of a non-bounding 1- 
cycle (a hole) that we refer to as a A-loop. From H10|) 
we can write down an integral for the probability, Pa (a) , 
that a PDC gives us a A-loop in the alpha-shape. Then, 
in a region R, the expected total number of holes at any 
radius a, is bounded by 

E0i(a) > A 2 ||P||P A (a), 

where A2 — 2A is the intensity of triangles in a 2D 
Poisson-Delaunay complex. 

The conditions on the size and shape of a triangle to 
generate A-loop are that 

1. all edges belong to the alpha-shape, i.e. Z max < 2a; 

2. the 2-simplex is excluded from the alpha-shape, i.e. 
the circumradius satisfies r > a; 

3. the circumcenter must be interior to the triangle, 
i.e. it is an acute triangle and the largest vertex 
angle satisfies max < 7r/2. 

The length of an edge in a triangle is related to the 
angle at the opposite vertex via I — 2rsin0. Thus, the 
condition Z max < 2a implies that r < a/ sin(0 max ). The 
marginal density for the largest angle at a vertex of a 
Poisson-Delaunay triangle is known to be [2^, p. 398] 



§ [(30 -tt) sin 20 

/max(0) = { - COS 20 + COS 40] , 

^ [sin + (tt — 0) cos 0] sin 0, 
Thus an integral expression for Pa (a) is 

fa/ sin (j) 



f < 



f < 



< - 

^ 2 



< TT. 



Pa = 



2( 7 rA)Ve- 7rAr / max (0)drd0. (12) 



To evaluate this integral, we start with an expression 
for the indefinite integral of the circumradius pdf: 



G(r) = / 2(ttA)V 



2^3„-7rAr- 



dr = -(n\r 2 + 1) 



-7rAr 



Evaluating with the limits of integration from l|12() we 
obtain 



G(t],<j>) = fa + l)e- 



V 



sin 



where we have simplified notation by using the reduced 
density, r\ = ir\a 2 . The second integral with respect to 
the angle does not have an analytic solution. However, 
we can obtain an approximate expression for small a (i.e. 
small rf) by using a Taylor expansion. First note that 

2 3 4 

/ t , X X X 

(* + i)e K = 1 -y + y-y + --- 

The first term in the series for G(rj, 0), and consequently 
the leading order term of Pa, is therefore rj 2 . The coef- 
ficient of rf , for j > 2 in the series for Pa is therefore 
given by the integral 



P 



(f) 



A 



1 



fn 



To evaluate these integrals requires only standard tech- 
niques from real calculus; for the first few terms we have: 

P A = 0.03200?/ 2 - 0.03422?7 3 + 0.01835?7 4 + 0{rf). 

Since the intensity of Delaunay cells is 2A, we have that 
the expectation per unit area of the first Betti number 
for small rj is 

E0i(v) > 2Xp A - 0.0640A?7 2 . 

In fact this lower bound on E0i(rj) is an asymptotic ex- 
pression as rj — > 0. This is because a connected cluster 
of at least three discs is needed to create a A-loop, and 
at least four overlapping discs are necessary to create a 
non-bounding 1-cycle with four or more edges. We know 
from the cluster expansions in Table [IJ however, that the 
leading order term as r/ — > for the number of fc-mers 
is r/ k ~ 1 . Thus there can be no other contribution to the 
rj 2 coefficient in a series expansion of Ef3\{rj). Indeed, in 
the limit of small 77, our simulations show exactly this 
behaviour — see Fig. 



C. Empty triangles in 3D 

We can derive a similar integral expression to that 
above for the probability of a A-loop in R 3 . However, 
in three dimensions not every A-loop represents an inde- 
pendent 1-cycle in the homology group. To see why this 
is the case, consider a cage consisting of the six edges of 
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a tetrahedron. There are four A-loops in this cage but 
only three independent 1-cycles, since the fourth A-loop 
is the sum of the other three. Nevertheless, using a simi- 
lar argument to that in the previous section, in the limit 
of small a, or small ij — ^nXa 3 , we can assume that the 
A-loops are isolated and that in a region R, 

Efoia) ~ A 2 ||i?||P A (a) 

where A 2 = ||7r 2 A is the intensity of faces in a 3D 
Poisson-Delaunay complex. 

The conditions for the existence of a A-loop in a 3D 
Poisson-Delaunay complex are essentially the same as 
those in two dimensions, except that they now apply to 
a typical face of a 3D PDC: 

1. all edges of the typical face belong to the alpha- 
shape, i.e. ^ max < 2a; 

2. the circumradius of the face satisfies p > a; 

3. the circumcenter of the face must be in the relative 
interior of the triangle, i.e., the largest vertex angle 
in a typical face satisfies max < ir/2. 

We use the relationship between edge-length and op- 
posite angle again so that condition 1 above becomes 
p < a I sin((/) max ). We also use the relationship between 
the face circumradius (p) and tetrahedron circumradius 
(r) of p = rsin^. Thus, an integral expression for P/\(a) 
in the three-dimensional setting is: 

a/ sin 4) p7r 



PA = I I / /maxWOMs^) SOT/eW^P#- 

Jf Ja JO 

The densities are 

1f)5 fl+c°s9 inc 

f e (0) = —j hsin 5 9dh = — sin 5 9(1 + cosfl) 2 , 



Where Z = (§7rAp 3 )/2, and M t ,M 2 are the Meijer G- 
functions: 



(4:1 * 



M 2 (z) = G$ [ z 



1 . _ 1 1 

2 ' 3 ' 3 

1 —I I I- 

2 ' 6 ' 6 ' 2 ' 



1.12 

2 ' 3 ' 3 

1 _I I 1- 

2 ' 6 ' 6 ' 2 ' ' 



Meijer G-functions are defined by integrals of Gamma 
functions 25 1. The form used within Mathematica is 




■A ) 



n™ x v(b 3 + s )m =1 v(i 



s) 



C L1 j= n +\ 



j=m+ 



X Y{\ - b 3 - s) 



ds. 



where the contour C divides the complex plane into two 
unbounded regions and separates the poles of r(l — a^ — s) 
and the poles of r(6, + s). 

Both M\(z) and M 2 (z) diverge as z — > 0. The products 
Z 3 Mi(Z 2 ) — > as Z — > 0, however, so we determine Tay- 
lor expansions about Z = for these terms. The zeroth 
and first order terms vanish and the second derivative 
has the value 



dZ 2 



[QZ^M^Z 2 ) + Z 3 M 2 (Z 2 



1 z=o 



Thus, to second order in 77 = |7rAa 3 



-35 A fri 



128V3 2 



(I)' 



13.8564 = A. 



1 



sin 



We can now compute the integral with respect to <p of 
H( r li4')fma.x(<t > ) and obtain the small r\ limit of 



fn 



sin 9 



sin 



exp 



ttA- 



sm 



^sin 2 </>[(3</>-7r)(3-2sin 2 

.. j — (9 — 16 sin 4 (/)) sin (/i cos < 

•''"- :U " - \ ^sin 2 0[(^-0)(3-2sin 2 ( / 

+ 3 cos 6 sin d>] 



■ e - - 

c 1-3 ' 2' 



The expression for / m ax(0) is due to Muche [22, p. 399]. 
We begin with the 9 integral: 



d9. 



F{P)= r-^(l + cos0) 2 expf-|7rA-^ 
Jo sm 9 \ sin 

An expression for F(p) may be given (using Mathemat- 
ica) in terms of Meijer G-functions and these are then 
integrated with respect to p to find 



01 r pa/sincf) 

Hic^^—HnXf F(p)dp 



-35 



128V3 



\§Z 3 M 1 (Z 2 ) + Z 3 M 2 (Z 



2\l a / s 



tt/2 



tt/3 



1 



sm 



-35 Arf 
128\/3 8 

35 A I 
12871 8 \^~4 ]V ~ 



- 1 /max (</>)<# 



The expectation per unit volume of the first Betti number 
for small rj is therefore 



35 



\/3A 
64 



7T 

T 



Ar/ 2 = 0.5747Ar/ 2 



This coefficient is exactly that obtained by comparing the 
expansions for \ and 0o, and agrees well with the value 
obtained in simulations. 



D. Empty tetrahedra in 3D 

Finally, we consider the existence of a 2-cycle in an 
alpha-shape formed by the four faces of a single Poisson- 
Delaunay tetrahedron. The conditions for this to occur 
are that 
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1. all faces belong to the alpha-shape, i.e. pi < a, 
where pi is the circumradius of the face opposite 
vertex Xf, 

2. the circumccntre is not covered by the union of balls 
of radius a, i.e. the tetrahedron circumradius sat- 
isfies r > a; 

3. the circumcenter must be interior to the tetrahe- 
dron. 

The conditions 1 and 3 both relate to the angle 9i, be- 
tween the outward-pointing normal to a face and a vector 
from the circumcenter to a vertex on that face. The cir- 
cumradius of face-i is pi — r sin 6^ , so condition 1 becomes 
r < a I sin 0j. The condition for the circumcenter to be 
interior to the tetrahedron requires that 9i < ir/2 for 
i = 0, 1,2,3. Thus an integral expression for an empty 
tetrahedron is: 



Ptet — 



7v/2 pa/ sin 9 



f m ,Mf R (r)drde, 



where fn{r) is the marginal pdf for the tetrahedron cir- 
cumradius defined in i|ll|) and f max (9) is an unknown pdf 
for the largest face-normal-vertex angle of a PDC. The 
lower limit, 9q, is the angle for a regular tetrahedron and 
is 9 — arccos | = 70.53°. 

Without knowing f ma , x (9), we can still find the leading 
order term for Ptet in the limit of small a. Firstly, the 
indefinite integral for the circumradius is 



G(r) 



^(§^A) 3 r 8 exp(-f^Ar 3 )dr 



+- ^TrAr 3 + 1 ] exp(-§7rAr 3 ) 




+ x + 1 ) exp(— x). 



FIG. 9: Blue dots show the normalised histogram for the 
typical face-normal - vertex angle obtained from a simulation 
of over four million tetrahedra. The data agree well with the 
known pdf for this quantity which is plotted as the pale-blue 
solid curve. The red dots mark the normalised histogram 
for the largest face-normal - vertex angle. We use this as a 
numerical approximation to /max(#). 



complex for a large number of points in a cube. The 
distribution of tetrahedra in a very large complex is ap- 
proximately the same as that obtained from many inde- 
pendent realisations. We generated 10 6 points with uni- 
form random coordinates in [—1, 1] , built the Delaunay 
complex, and discarded tetrahedra with circumcenters 
within a 0.2 margin of the boundary to minimize edge 
effects. This yielded over four million Poisson-Delaunay 
tetrahedra. The probability density for the typical face- 
normal-vertex angle is known to be |23| 



where we have simplified notation by using the reduced 
density, x = |7rAr 3 . The Taylor expansion for small x is 



-{\x 2 + x + l)e x 



x 

+ ¥ 



0(x 4 ), 



so that to highest order in r/ = ^nXa 3 , 
G(r h 9) 



rr_ 
6 



In terms of Pt e t we have 

„3 fi/2 



Pi 



lei 



1 dO. 



(13) 



We are unable to derive an analytic expression for 
/max(#), so we estimate it by simulation and calculate 
a numerical approximation to Ptet- The Poisson point 
process is ergodic, so the simplest technique for simu- 
lating Poisson-Delaunay cells is to build the Delaunay 



1 0S 

f(0) = — sin 5 0(1 



COS! 



and provides a check on our simulation. Normalised 
histograms for the typical and the largest face-normal- 
vertex angle in a PDC are shown in Fig. [!|] The numer- 
ical approximation to the integrand in (|13|l is shown in 
Fig. ^3 The area under the curve as calculated from this 
data is 0.0023. Now, since the intensity of Delaunay cells 
is A3 = ||7r 2 A, we have that the expectation per unit 
volume of the second Betti number in the limit of small 
r\ is 



24 

Efo(v) - ^ 2 ^ p tet - 0.015 A77 3 . 
35 

Again, we know this result is an asymptotic one because 
2-cycles that involve more than the faces of a single De- 
launay tetrahedron are necessarily built from five or more 
overlapping balls. As we see in Table [HI the expected 
number of such clusters has leading order 77 as rj — > 0. 
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0.014 F ' 1 1 1 1 ' = 

0.012- ■■>'' 
0.01 - 

0.008 - '; *. 

0.006 - ■■• \ 
0.004 - 
0.002 - 

o ~ ,nr ''" 1 1 1 ' 1 

1.25 1.3 1.35 1.4 1.45 1.5 1.55 

FIG. 10: Blue dots show the numerical approximation to 
the integrand in Eq. 11311 over the domain of integration 
arccos § < < f • 



E. Further analysis 

In the study of percolation theory or coverage pro- 
cesses, the total number of connected components is stud- 
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